Wind trajectory
library(GeoPressureR)
library(tidyverse)
library(leaflet)
library(leaflet.extras)
library(raster)
library(dplyr)
library(ggplot2)
library(plotly)
knitr::opts_chunk$set(echo = FALSE)
load(paste0("../data/1_pressure//", params$gdl_id, "_pressure_prob.Rdata"))
load(paste0("../data/3_static/", params$gdl_id, "_static_prob.Rdata"))
# load(paste0("../data/4_basic_graph/", params$gdl_id, "_basic_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_wind_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_grl.Rdata"))
col <- rep(RColorBrewer::brewer.pal(8, "Dark2"), times = ceiling(max(pam$sta$sta_id) / 8))
Altitudes are computed based on pressure measurement of the geolocation, corrected based on the assumed location of the shortest path. This correction accounts therefore for the natural variation of pressure as estimated by ERA-5. The vertical lines indicate the sunrise (dashed) and sunset (solid).
p <- ggplot() +
# geom_line(data = pam$pressure, aes(x = date, y = obs), colour = "grey") +
geom_line(data = do.call("rbind", shortest_path_timeserie), aes(x = date, y = altitude)) +
geom_line(data = do.call("rbind", shortest_path_timeserie) %>% filter(sta_id > 0), aes(x = date, y = altitude, col = factor(sta_id))) +
# geom_vline(data = twl, aes(xintercept = twilight, linetype = ifelse(rise, "dashed", "solid"), color="grey"), lwd=0.1) +
theme_bw() +
scale_colour_manual(values = col) +
scale_y_continuous(name = "Altitude (m)")
ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)
file <- paste0("figure_print/wintering_location/wintering_location_",params$gdl_id,".png")
if(file.exists(file)){
knitr::include_graphics(file)
}
tmp <- lapply(pressure_prob, function(x) {
mt <- metadata(x)
df <- data.frame(
start = mt$temporal_extent[1],
end = mt$temporal_extent[2],
sta_id = mt$sta_id
)
})
tmp2 <- do.call("rbind", tmp)
sim_lat <- as.data.frame(t(path_sim$lat)) %>%
mutate(sta_id = path_sim$sta_id) %>%
pivot_longer(-c(sta_id)) %>%
left_join(tmp2,by="sta_id")
sim_lat_p <- sim_lat %>%
filter(sta_id==max(sta_id)) %>%
mutate(start=end) %>%
rbind(sim_lat)
sp_lat <- as.data.frame(shortest_path) %>% left_join(tmp2,by="sta_id")
sp_lat_p <- sp_lat %>%
filter(sta_id==max(sta_id)) %>%
mutate(start=end) %>%
rbind(sp_lat)
p <- ggplot() +
geom_step(data=sim_lat_p, aes(x=start, y=value, group=name), alpha=.07) +
geom_point(data=sp_lat_p, aes(x=start, y=lat)) +
xlab('Date') +
ylab('Latitude') +
theme_light()
ggplotly(p, dynamicTicks = T)
The large circles indicates the shortest path (overall most likely trajectory) estimated by the graph approach. The size is proportional to the duration of stay. The small dots and grey lines represents 10 possible trajeectories of the bird according to the model.
Click on the full-screen mode button on the top-left of the map to see more details on the map.
sta_duration <- unlist(lapply(static_prob_marginal, function(x) {
as.numeric(difftime(metadata(x)$temporal_extent[2], metadata(x)$temporal_extent[1], units = "days"))
}))
pal <- colorFactor(col, as.factor(seq_len(length(col))))
m <- leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>%
addFullscreenControl() %>%
addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = pal(factor(shortest_path$sta_id, levels = pam$sta$sta_id)), weight = sta_duration^(0.3) * 10)
for (i in seq_len(nrow(path_sim$lon))) {
m <- m %>%
addPolylines(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = 0.5, weight = 1, color = "#808080") %>%
addCircles(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = .7, weight = 1, color = pal(factor(shortest_path$sta_id, levels = pam$sta$sta_id)))
}
m
The marginal probability map estimate the overall probability of position at each stationary period regardless of the trajectory taken by the bird. It is the most useful quantification of the uncertainty of the position of the bird.
li_s <- list()
l <- leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>%
addFullscreenControl()
for (i_r in seq_len(length(static_prob_marginal))) {
i_s <- metadata(static_prob_marginal[[i_r]])$sta_id
info <- metadata(static_prob_marginal[[i_r]])$temporal_extent
info_str <- paste0(i_s, " | ", info[1], "->", info[2])
li_s <- append(li_s, info_str)
l <- l %>%
addRasterImage(static_prob_marginal[[i_r]], colors = "OrRd", opacity = 0.8, group = info_str) %>%
addCircles(lng = shortest_path$lon[i_s], lat = shortest_path$lat[i_s], opacity = 1, color = "#000", weight = 10, group = info_str)
}
l %>%
addLayersControl(
overlayGroups = li_s,
options = layersControlOptions(collapsed = FALSE)
) %>%
hideGroup(tail(li_s, length(li_s) - 1))
fun_marker_color <- function(norm){
if (norm < 20){
"darkpurple"
} else if (norm < 35){
"darkblue"
} else if (norm < 50){
"lightblue"
} else if (norm < 60){
"lightgreen"
} else if (norm < 80){
"yellow"
} else if (norm < 100){
"lightred"
} else {
"darkred"
}
}
fun_NSEW <- function(angle){
angle <- angle %% (pi* 2)
angle <- angle*180/pi
if (angle < 45/2){
"E"
} else if (angle < 45*3/2){
"NE"
} else if (angle < 45*5/2){
"N"
} else if (angle < 45*7/2){
"NW"
} else if (angle < 45*9/2){
"W"
} else if (angle < 45*11/2){
"SW"
} else if (angle < 45*13/2){
"S"
}else if (angle < 45*15/2){
"SE"
} else {
"E"
}
}
sta_duration <- unlist(lapply(static_prob_marginal,function(x){as.numeric(difftime(metadata(x)$temporal_extent[2],metadata(x)$temporal_extent[1],units="days"))}))
m <-leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>% addFullscreenControl() %>%
addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#000", weight = sta_duration^(0.3)*10)
for (i_s in seq_len(grl$sz[3]-1)){
if (grl$flight_duration[i_s]>5){
edge <- which(grl$s == shortest_path$id[i_s] & grl$t == shortest_path$id[i_s+1])
label = paste0( i_s,': ', grl$flight[[i_s]]$start, " - ", grl$flight[[i_s]]$end, "<br>",
"F. dur.: ", round(grl$flight_duration[i_s]), ' h <br>',
"GS: ", round(abs(grl$gs[edge])), ' km/h, ',fun_NSEW(Arg(grl$gs[edge])),'<br>',
"WS: ", round(abs(grl$ws[edge])), ' km/h, ',fun_NSEW(Arg(grl$ws[edge])),'<br>',
"AS: ", round(abs(grl$as[edge])), ' km/h, ',fun_NSEW(Arg(grl$as[edge])),'<br>')
iconArrow <- makeAwesomeIcon(icon = "arrow-up",
library = "fa",
iconColor = "#FFF",
iconRotate = (90 - Arg(grl$ws[edge])/pi*180) %% 360,
squareMarker = TRUE,
markerColor = fun_marker_color(abs(grl$ws[edge])))
m <- m %>% addAwesomeMarkers(lng = (shortest_path$lon[i_s] + shortest_path$lon[i_s+1])/2,
lat = (shortest_path$lat[i_s] + shortest_path$lat[i_s+1])/2,
icon = iconArrow, popup = label)
}
}
m
edge <- t(graph_path2edge(path_sim$id, grl))
nj <- ncol(edge)
nsta <- ncol(path_sim$lon)
speed_df <- data.frame(
as = abs(grl$as[edge]),
gs = abs(grl$gs[edge]),
ws = abs(grl$ws[edge]),
sta_id_s = rep(head(grl$sta_id,-1), nj),
sta_id_t = rep(tail(grl$sta_id,-1), nj),
flight_duration = rep(head(grl$flight_duration,-1), nj),
dist = geosphere::distGeo(
cbind(as.vector(t(path_sim$lon[,1:nsta-1])), as.vector(t(path_sim$lat[,1:nsta-1]))),
cbind(as.vector(t(path_sim$lon[,2:nsta])), as.vector(t(path_sim$lat[,2:nsta])))
) / 1000
) %>% mutate(
name = paste(sta_id_s,sta_id_t, sep="-")
)
plot1 <- ggplot(speed_df, aes(reorder(name, sta_id_s), gs)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot2 <- ggplot(speed_df, aes(reorder(name, sta_id_s), ws)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot3 <- ggplot(speed_df, aes(reorder(name, sta_id_s), as)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot4 <- ggplot(speed_df, aes(reorder(name, sta_id_s), flight_duration)) + geom_point() + theme_bw() +scale_x_discrete(name = "")
subplot(ggplotly(plot1), ggplotly(plot2), ggplotly(plot3), ggplotly(plot4), nrows=4, titleY=T)
alt_df = do.call("rbind", shortest_path_timeserie) %>%
arrange(date) %>%
mutate(
sta_id_s = cummax(sta_id),
sta_id_t = sta_id_s+1
) %>%
filter(sta_id == 0 & sta_id_s > 0 ) %>%
group_by(sta_id_s, sta_id_t) %>%
summarise(
alt_min = min(altitude),
alt_max = max(altitude),
alt_mean = mean(altitude),
alt_med = median(altitude),
)
trans_df <- speed_df %>%
group_by(sta_id_s,sta_id_t,flight_duration) %>%
summarise(
as_m = mean(as),
as_s = sd(as),
gs_m = mean(gs),
gs_s = sd(gs),
ws_m = mean(ws),
ws_s = sd(ws),
dist_m = mean(dist),
dist_s = sd(dist)
) %>%
left_join(alt_df)
trans_df %>% kable()
| sta_id_s | sta_id_t | flight_duration | as_m | as_s | gs_m | gs_s | ws_m | ws_s | dist_m | dist_s | alt_min | alt_max | alt_mean | alt_med |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 3 | 25.0 | 34.78721 | 1.561995 | 48.89101 | 1.666275 | 25.06176 | 0.7228002 | 1222.27520 | 41.65688 | NA | NA | NA | NA |
| 3 | 4 | 1.0 | 36.24353 | 17.683192 | 42.10700 | 26.717529 | 27.45286 | 1.4918548 | 42.10700 | 26.71753 | 139.19808 | 277.4258 | 208.4051 | 208.5913 |
| 4 | 5 | 2.0 | 30.48163 | 12.644211 | 30.89662 | 13.183696 | 21.91755 | 1.8785415 | 61.79325 | 26.36739 | 101.58665 | 304.4136 | 175.5619 | 137.1704 |
| 5 | 6 | 1.5 | 42.75670 | 18.500498 | 32.64181 | 19.142397 | 22.15592 | 4.1416911 | 48.96272 | 28.71360 | 57.92701 | 396.5567 | 183.9772 | 140.7125 |
| 6 | 7 | 26.5 | 35.29000 | 5.796444 | 44.23211 | 3.363960 | 34.42406 | 2.5794491 | 1172.15101 | 89.14494 | -88.36091 | 1397.1578 | 308.1802 | 128.1707 |
| 7 | 8 | 6.0 | 55.01019 | 12.787171 | 81.05204 | 19.244547 | 62.06125 | 4.5192942 | 486.31222 | 115.46728 | 92.94742 | 1255.0715 | 508.4330 | 285.3609 |
| 8 | 9 | 7.0 | 42.02473 | 15.230921 | 30.30106 | 18.132733 | 20.17888 | 1.3850864 | 212.10742 | 126.92913 | 26.55118 | 1060.9777 | 218.3246 | 108.0140 |
| 9 | 10 | 3.0 | 34.07205 | 16.007661 | 40.73461 | 19.981061 | 16.45508 | 4.9016519 | 122.20382 | 59.94318 | NA | NA | NA | NA |
| 10 | 11 | 4.0 | 39.62759 | 17.378455 | 44.56247 | 22.262268 | 15.98501 | 3.8274602 | 178.24990 | 89.04907 | 325.47014 | 1092.1917 | 645.5386 | 580.9072 |
| 11 | 12 | 3.0 | 34.62132 | 15.678381 | 37.51698 | 14.527980 | 34.71207 | 1.6842107 | 112.55093 | 43.58394 | 115.51513 | 911.4002 | 473.4592 | 296.2463 |
| 12 | 13 | 2.5 | 31.67524 | 10.830975 | 27.94943 | 16.033425 | 28.34879 | 1.2195439 | 69.87357 | 40.08356 | 61.25517 | 877.9342 | 450.4868 | 417.1281 |